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1. Introduction 

In the last few years solid state physics has increasingly benefited from scientific 
computing and the significance of numerical techniques is likely to keep on growing 
quickly in this field. Because of the high complexity of solids, which are made of a huge 
number of interacting electrons and nuclei, a full understanding of their properties cannot 
be developed using analytical methods only. Numerical simulations do not only provide 
quantitative results for the properties of specific materials but are also widely used to 
test the validity of theories and analytical approaches. 

Numerical and analytical approaches based on perturbation theory and effective 
independent-particle theories such as the Fermi liquid theory, the density functional 
theory, the Hartree-Fock approximation, or the Born-Oppenheimer approximation, have 
been extremely successful in explaining the properties of solids. However, the low-energy 
and low-temperature electronic, optical, or magnetic properties of various novel ma- 
terials are not understood within these simplified theories. For example, in strongly 
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correlated systems, the interactions between constituents of the sohd are so strong that 
they can no longer be considered separately and a collective behavior can emerge. As 
a result, these systems may exhibit new and fascinating macroscopic properties such 
as high-temperature superconductivity or colossal magneto-resistance Quasi-one- 
dimensional electron-phonon (EP) systems like MX-chain compounds are other exam- 
ples of electronic systems that are very different from traditional ones [2]. Their study 
is particularly rewarding for a number of reasons. First they exhibit a remarkably wide 
range of competing forces, which gives rise to a rich variety of different phases charac- 
terized by symmetry-broken ground states and long-range orders. Second these systems 
share fundamental features with higher-dimensional novel materials (for instance, high- 
Tc cuprates or charge-ordered nickelates) such as the interplay of charge, spin, and lattice 
degrees of freedom. One-dimcnsional (ID) models allow us to investigate this complex 
interplay, which is important but poorly understood also in 2D- and 3D highly correlated 
electron systems, in a context more favorable to numerical simulations. 

I'l. Models. - Calculating the low-energy low-temperature properties of solids from 
first principles is possible only with various approximations which often are not reliable 
in strongly correlated or low-dimensional electronic systems. An alternative approach 
for investigating these materials is the study of simplified lattice models which include 
only the relevant degrees of freedom and interactions but nevertheless are believed to 
reproduce the essential physical properties of the full system. 

A fundamental model for ID correlated electronic systems is the Hubbard model |3] 
defined by the Hamiltonian 



It describes electrons with spin a =T,i which can hop between neighboring sites on a 
lattice. Here c|^, c^^ are creation and annihilation operators for electrons with spin a at 
site i, riicr — cl^c^^ ^^'^ ^^'^ corresponding density operators. The hopping integral t gives 
rise to a a single-electron band of width AtD (with D being the spatial dimension) . The 
Coulomb repulsion between electrons is mimicked by a local Hubbard interaction U > 0. 
The average electronic density per site is < n < 2, where n ~ N^/N, N is the number 
of lattice sites and is the number of electrons; n/2 is called the band filling. 

In the past the Hubbard model was intensively studied with respect to (itinerant) 
ferromagnetism, antiferromagnetism and metal-insulator (Mott) transitions in transition 
metals. More recently it has been used in the context of heavy fermions and high- 
temperature superconductivity as perhaps the most fundamental model accounting for 
strong electronic correlation effects in solids. 

The coupling between electrons and the lattice relaxation and vibrations (phonons) 
is also known to have significant effects on the properties of solids including the above 
mentioned strongly correlated electronic materials (see the papers by Egami, Calvani, 
Zhou, Perroni, and Saini). Dynamical phonon effects are particularly important in quasi- 
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Fig. 1. - Schematic representation of the ID Holstein- Hubbard model. 



ID metals and chargc-dcnsity-wavc (CDW) systems. The simplest model describing the 
effect of an additional EP coupling is the Holstein-Hubbard model. This model describes 
electrons coupled to dispersionless phonons, represented by local Einstein oscillators. 
The Hamiltonian is given by 

(2) -f^ep = ii^oo + ^Xl^^ + tXI^? ~ "H*"*' 

i i i 

where qi and pi are the position and momentum operators for a phonon mode at site i, 
and rii = + n^j^. At first sight, there are three additional parameters in this model 
(compared to the purely electronic model): The oscillator mass M, the spring constant 
A', and the EP coupling constant a. However, if we introduce phonon (boson) creation 
and annihilation operators b\ and b^ , respectively, the Holstein-Hubbard Hamiltonian 
can be written (up to a constant term) 

(3) ifep = i?cc + ^0 ^ b\b^ - gujo (4 + h K , 

i i 

where the phonon frequency is given by Uq = K/M (fi = 1) and a dimensionless EP 
coupling constant is defined hy g = aajujQ with the range of zero-point fluctuations 
given by 2a^ ~ (XA/)~^/^. We can set the parameter a equal to 1 by redefining the 
units of oscillator displacements. Thus, the effects of the EP coupling are determined 
by two dimensionless parameter ratios only: Lo^jt and g. Alternatively, one can use the 
polaron binding energy ~ g^i^o or, equivalently, A = ep/2tD, instead of g. The various 
constituents and couplings of the Holstein-Hubbard model are summarized in fig. ^ 

In the single-electron case, the resulting Holstein model 0] has been studied as a 
paradigmatic model for polaron formation (see the paper by Fehske, Alvermann, Ho- 
henadler and Wellein). At half-filling, the EP coupling may lead to a Peierls instability 
related to the appearance of CDW order in competition with the SDW instability trig- 
gered by U (see the separate paper by Fehske and Jeckelmann). 
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1'2. Methods. - Despite the great simplification brought by the above models, the- 
oretical investigations remain difRcult because a quantum many-particle problem has 
to be solved with high accuracy to determine correlation effects on physical properties 
beyond the mean-field level. Analytical solutions of these models are known for special 
cases only. To determine the spectral and thermodynamical properties of these models, 
theorists have turned to numerical simulations. Among the various approaches, exact 
numerical methods play a significant role. A numerical calculation is said to be exact 
if no approximation is involved aside from the restriction imposed by finite computa- 
tional resources (in particular, the same calculation would be mathematically exact if 
it were carried out analytically), the accuracy can be systematically improved with in- 
creasing computational effort, and actual numerical errors are quantifiable or completely 
negligible. Especially for strongly correlated systems, exact numerical methods are of- 
ten the only approach available to obtain accurate quantitative results in model systems 
and the results that they provide are essential for checking the validity of theories or 
testing the accuracy of approximative analytical methods. Nowadays, finite-cluster ex- 
act diagonalizations (ED), the numerical renormalization group (NRG), density matrix 
renormalization group (DMRG) calculations, quantum Monte Carlo (QMC) simulations, 
and the dynamical mean-field theory (DMFT) have become very powerful and impor- 
tant tools for solving many-body problems. In what follows, we briefly summarize the 
advantages and weaknesses of these techniques, especially in the context of EP lattice 
models such as the the Holstein-Hubbard model. Then in the remaining sections we will 
present the basic principles of the ED (sect.|2l) and DMRG (sect. |3| and |4|) approaches. 

The paradigm of an exact numerical calculation in a quantum system is the exact 
diagonalization (ED) of its Hamiltonian, which can be carried out using several well- 
established algorithms (sec the next section). ED techniques can be used to calculate 
most properties of any quantum system. Unfortunately, they are restricted to small 
systems (for instance, 16 sites for a half-filled Hubbard model and less for the Holstein- 
Hubbard model) because of the exponential increase of the computational effort with 
the number of particles. In many cases, these system sizes are too small to simulate 
adequately macroscopic properties of solids (such as (band) transport, magnetism or 
charge long-range order). They are of great value, however, for a description of local or 
short-range order effects. 

In QMC simulations the problem of computing quantum properties is transformed 
into the summation of a huge number of classical variables, which is carried out using 
statistical (Monte Carlo) techniques. There are numerous variations of this principle also 
for EP problems (see, e.g., the paper on QMC by Mishchenko). QMC simulations are 
almost as widely applicable as ED techniques but can be applied to much larger systems. 
The results of QMC calculations are affected by statistical errors which, in principle, 
can be systematically reduced with increasing computational effort. Therefore, QMC 
techniques arc often numerically exact. In practice, several problems such as the sign 
problem (typically, in frustrated quantum systems) or the large auto-correlation time in 
the Markov chain (for instance, in critical systems) severely limit the applicability and ac- 
curacy of QMC simulations in strongly correlated or low-dimensional systems. Moreover, 
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real-frequency dynamical properties often have to be calculated from the imaginary-time 
correlation functions obtained with QMC using an analytic continuation. However, the 
transformation of imaginary-time data affected by statistical errors to real-frequency 
data is an ill-conditioned numerical problem. In practice, one has to rely on approxi- 
mate transformations using least-square or maximum-entropy fits, which yield results of 
unknown accuracy. 

DMRG methods allow us to calculate the static and dynamical properties of quantum 
systems much larger than those possible with ED techniques (see the third and fourth 
section, as well as the separate paper by Fchske and Jeckelmann). For ID systems and 
quantum impurity problems, one can simulate lattice sizes large enough to determine 
static properties in the thermodynamic limit and the dynamical spectra of macroscopic 
systems exactly. In higher dimensions DMRG results are numerically exact only for sys- 
tem sizes barely larger than those available with ED techniques. For larger system sizes 
in dimension two and higher DMRG usually provides only a variational approximation 
of the system properties. 

The NRG is a precursor of the DMRG method. Therefore, it it not surprising that 
most NRG calculations can also be carried out with DMRG. NRG provides numerically 
exact results for the low-energy properties of quantum impurity problems (see the paper 
on NRG methods by Hewson) . Moreover, for this type of problem it is computationally 
more efficient than DMRG. For other problems (lattice problems, high-energy properties) 
NRG usually fails or provides results of poor accuracy. 

In the dynamical mean-field theory (see the papers by Ciuchi, Capone, and Castel- 
lani) it is assumed that the self-energy of the quantum many-body system is momentum- 
independent. While this is exact on a lattice with an infinitely large coordination number 
(i.e., in the limit of infinite dimension), it is considered to be a reasonable approxima- 
tion for 3D systems. Thus in applications to real materials DMFT is never an exact 
numerical approach, although a DMFT-based approach for a 3D system could possibly 
yield better results than direct QMC or DMRG calculations. It should be noticed that in 
the DMFT framework the self-energy has to be determined self-consistently by solving a 
quantum impurity problem, which is itself a difficult strongly correlated problem. This 
quantum impurity problem is usually solved numerically using one of the standard meth- 
ods discussed here (ED, NRG, QMC or DMRG). Therefore, the DMFT approach and 
its extensions can be viewed as an (approximate) way of circumventing the limitations 
of the standard methods and extend their applicability to large 3D systems. 

In summary, every numerical method has advantages and weaknesses. ED is unbiased 
but is restricted to small clusters. For ID systems DMRG is usually the best method while 
for 3D systems only direct QMC simulations are possible. QMC techniques represent 
also the most successful approach for non-critical and non-frustrated two-dimensional 
systems. There is currently no satisfactory numerical (or analytical) method for two- 
dimensional strongly correlated systems with frustration or in a critical regime. 
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2. Exact diagonalization approach 

As stated above, ED is presently probably the best controlled numerical method be- 
cause it allows an approximation-free treatment of coupled electron-phonon models in the 
whole parameter range. As a precondition we have to work with finite systems and apply 
a well-defined truncation procedure for the phonon sector (see subsect. l2"l|l . At least for 
the single-electron Holstein model a variational basis can be constructed in such a way 
that the ground-state properties of the model can be computed numerically exact in the 
thermodynamic limit (cf. subscct. 12 '2(1 . In both cases the resulting numerical problem 
is to find the eigenstates of a (sparse) Hermitian matrix using Lanczos or other iterative 
subspace methods (subsect. 12 '3f) . In general the computational requirements of these 
eigenvalue algorithms are determined by matrix- vector multiplications (MVM), which 
have to be implemented in a parallel, fast and memory saving way on modern super- 
computers. Extensions for the calculation of dynamical quantities have been developed 
on the basis of Lanczos recursion and kernel polynomial expansions (cf. subsect. I2'4|l . 
Quite recently cluster perturbation theory (CPT) has been used in combination with 
these techniques to determine the single-particle electron and phonon spectra. 

2'1. Many-body Hilbert space and basis construction. - 

2 '1.1. Basis symmetrization. The total Hilbert space of Holstein-Hubbard type mod- 
els (01 can be written as the tensorial product space of electrons and phonons, spanned 
by the complete basis set {\b) = |e) ® \p)} with 

N N 

(4) |e)-n n (4r-|0)e and |p) = [] ^=(6t)"v.|o)^. 

Here n^o-.e G {0, 1}, i.e. the electronic Wannier site i might be empty, singly or doubly oc- 
cupied, whereas we have no such restriction for the phonon number, nii p g {0, . . . , oo}. 
Consequently, e = 1, . . . ,De and p = 1, . . . ,Dp label basic states of the electronic and 
phononic subspaces having dimensions De = [J^ ){n^ ) ~ I'cspectively. 

Since the Holstein Hubbard Hamiltonian commutes with the total electron number op- 
erator Ne = X]t=i("i,T + "-i-i); ^e,cr = '^iLi^i.cr (we used the 'hat' to discriminate 
operators from the corresponding particle numbers), and the z-component of the total 
spin 5*^ = i X]i!=i('^i,T ~ basis {\b)} has been constructed for fixed and S^. 

To further reduce the dimension of the total Hilbert space, we can exploit the space 
group symmetries [translations (Gt) and point group operations (Gl)] and the spin-flip 
invariance [(Gs); = - subspace only]. Clearly, working on finite bipartite clusters 
in ID or 2D (here N = k'^ + P, k and I are both even or odd integers) with periodic 
boundary conditions (PBC) , we do not have all the symmetry properties of the underlying 
ID or 2D (square) lattices Restricting ourselves to the ID non-cquivalcnt irreducible 
representations of the group G{K) = Gt x Gl{K) x Gs, wc can use the projection 
operator T'^^, = [.9(^)1 EeeG(i?) xg,^, Q (with [H^V^ J = 0, V^^^^ = V^ ^^ and 
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Vj^ rs ^i?' r's' ~ '^K rs K' ^^.r' ^s,s') in Order to generate a new symmetrized basis set: 

{\b)} —* {|6)}. Q denotes the g{K) elements of the group G{K) and X'^\^ is the (complex) 

character of Q in the [K, rs] representation, where K refers to one of the N allowed wave 
vectors in the first Brillouin zone, r labels the irreducible representations of the little 
group of K, Gl{K), and s parameterizes Gs- For an efficient parallel implementation 
of the MVM it is extremely important that the symmetrized basis can be constructed 
preserving the tensor product structure of the Hilbert space, i.e., 

(5) {\b)=Nf'-'^rg^^, [\e)®\p)]}, 



with e = 1, . . . , De [De ^ De/g{K)]. The A'^ arc normalization factors. 

2'1.2. Phonon Hilbert space truncation. Since the Hilbert space associated to the 
phonons is infinite even for a finite system, we apply a truncation procedure [HI retaining 
only basis states with at most M phonons: 

N 

(6) {\p) ; nip = ^m,,p < M}. 

i=l 

The resulting Hilbert space has a total dimension D = De^^^ x D^^ with D^^ = ^^j^^r', 
and a general state of the Holstein Hubbard model is represented as 

e— 1 p—1 

It is worthwhile to point out that, switching from a real-space representation to a momen- 
tum space description, our truncation scheme takes into account all dynamical phonon 
modes. This has to be contrasted with the frequently used single- mode approach 0. In 
other words, depending on the model parameters and the band filling, the system "de- 
cides" by itself how the M phonons will be distributed among the independent Einstein 
oscillators related to the N Wannier sites or, alternatively, among the N different phonon 
modes in momentum space. Hence with the same accuracy phonon dynamical effects on 
lattice distortions being quasi-localized in real space (such as polarons, Frenkel excitons, 
. . . ) or in momentum space (like charge-density- waves, . . . ) can be studied. 

Of course, one has carefully to check for the convergence of the above truncation 
procedure by calculating the ground-state energy as a function of the cut-off parameter 
M . In the numerical work convergence is assumed to be achieved if Eq is determined 
with a relative error A£;^^^^ = {Eo{M) - Eq{M - 1))/Eq{M) < 10"^ . In addition we 
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Fig. 2. - Convergence of the phonon distribution function |c'"p(M) and ground-state energy 
Eo{M) (inset; here filled symbols give the results obtained separating the Q = phonon mode) 
as a function of the maximum number of phonons M retained. Results are given for the Holstein 
model on a ID lattice with N = 14 (a) and A'^ = 6 (b) sites (PBC), where the parameters Sp 
and ujQ are the same as for the CPT calculation presented in fig. 3 of the paper by Fehske, 
Alvermann, Hohenadler and Wellein. 



guarantee that the phonon distribution function 

(8) |cMp(Af)= J2 E Ic^pI'' 

e=l p=l 

{mp=m} 

which gives the different weights of the m-phonon states in the ground-state IV'o) i becomes 
independent of AI and |c'^^^p(M) < 10^®. To illustrate the M dependences of the 
phonon distribution function and the ground-state energy, we have shown both quantities 
in fig. 121 for the single-electron Holstein model on rather small lattices. Figure |21 proves 
that our truncation procedure is well controlled even in the strong EP coupling regime, 
where multi-phonon states become increasingly important. 

For the Holstein-type models the computational requirements can be further reduced. 
Here it is possible to separate the symmetric phonon mode, Bq = bi, and to calcu- 

late its contribution to H analytically jS] . For the sake of simplicity, we restrict ourselves 
to the ID spinless case in what follows. Using the momentum space representation of 
the phonon operators, the original Holstein Hamiltonian takes the form 

(9) H = -tY^icU, + c]q) - E(^-Q. + ^Q>Q. + E ^L^Q. 

ij 3 3 

^ith B]^^ = Uj^ih\, Bq^ = Ulibi = and uq^ = J2^^3,^r>'^^ where U^,, = (l/\/iV)x 

expliQjRi} and Qj (Ri) denote the allowed momentum (translation) vectors of the 
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Fig. 3. - Small variational Hilbert space for the 
ID polaron problem. Basis states are repre- 
sented by dots, off-diagonal matrix elements by 
lines. Vertical bonds create or destroy phonons 
with frequency uq. Horizontal bonds correspond 
to electron hops (oc t). Accordingly, state |1) 
describes an electron at the origin (0) and no 
phonon, state |2) is an electron and one phonon 
both at site 0, |3) is an electron at the nearest 
neighbor-site site 1, and a phonon at site 0, and 
so on. The figure is re-drawn from ref. 



lattice. The Q = phonon mode couples to tlq ~ Nf,/ \/N which is a constant if working 
in a subspace with fixed number of electrons. Thus the Hamiltonian decomposes into 
H = H'+Hq^q, with Hq^o = -^£^ {bI+B„) Uo+ujoBIBo. Since [H' , Hq=o] = 0, the 
eigenspectrum of H can be built up by the analytic solution for Hq^q and the numerical 
results for H' . Using the unitary transformation 

(10) Sm = e.,[-^^^iBl-B,)] 

which introduces a shift of the phonon operators {Bq ^ Bq + ^J^) , we easily find the 
diagonal form of Hq=o ~ ojqBqBq — SpN'^/N. It represents a harmonic oscillator with 
eigenvalues and eigenvectors Ej ~ luqI — and \l) ~ -^(_bJ)'^|0). The corresponding 

eigenenergies and eigenvectors of Hq^q are Ei = Ej and \l{Ne)) = S^{Ne)\l), respectively. 
That is, in the eigenstates of the Holstein model a homogeneous lattice distortion occurs. 
Note that the homogeneous lattice distortions are different in subspaces with different 
electron number. Thus excitations due to lattice relaxation processes will show up in 
the one-particle spectral function. Finally, eigenvectors and eigenenergies of H can be 
constructed by combining the above analytical result with the numerically determined 
eigensystem {E',^; \^'^)) of H': E'^ ^ = ^ + uo^l - EpN^/N and \i'n,i) = |0 «> l^(^e)). 

2'2. Variational ED method. ~ In this section we briefly outline a very efficient varia- 
tional method to address the (one-electron) Holstein polaron problem numerically in any 
dimension. The approach was developed by Bonca et al. |51 11U| and is based on a clever 
way of constructing the EP Hilbert space which can be systematically expanded in order 
to achieve high-accuracy results with rather modest computational resources. 

The authors built up the variational space starting from an initial state, e.g. the 
electron at the origin, and acting repeatedly {L times) with the off-diagonal diagonal 
hopping [t] and EP coupling (A) terms of the Hamiltonian (see fig. O. A basis state 
is added if it is connected by a non-zero t- or A-matrix element to a state previously 
in the space, i.e., states in generation I are obtained by acting / times with off-diagonal 
terms. Only one copy of each state is retained. Importantly, all translations of these 
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states on an infinite lattice are included. According to Bloch's theorem each eigenstate 
can be written a,s ip = b^'^^ql, where is a set of complex amplitudes related to the 
states in the unit cell, e.g. L = 7 for the small variational space shown in fig. For 
each momentum K the resulting numerical problem is then to diagonalize a Hermitian 
L X L matrix. While the size of the Hilbert space increases as (D+1)-'^, the error in the 
ground-state energy decreases exponentially with L. Thus in most cases 10'*-10® basis 
states are sufficient to obtain an 8-10 digit accuracy for Eq. The ground-state energy 
calculated this way is variational for the infinite system. 

2'3. Solving the eigenvalue problem. - To determine the eigenvalues of large sparse 
Hermitian matrices, iterative (Krylov) subspace methods like Lanczos |11| and variants 
of Davidson |12| diagonalization techniques are frequently applied. These algorithms 
contain basically three steps: 

(1) project problem matrix A £ R" onto a subspace A*^ G V*"' (fc <C n) 

(2) solve the eigenvalue problem in V'^ using standard routines 

(3) extend subspace y/^+i by ^ vector t _L V*^ and go back to (2). 

This way we obtain a sequence of approximative inverses of the original matrix A. 

2'3.1. Lanczos diagonalization. Starting out from an arbitrary (random) initial state 
|(/9o)7 having finite overlap with the true ground state lipo): the Lanczos algorithm recur- 
sively generates a set of orthogonal states (Lanczos or Krylov vectors) : 

(11) =n^\^i)-ai\^i)-b^\^i^,), 

where a/ = {(fi\lI°\Lpi)/{ipi\(fi),bf = ((^/|(^;)/(<^;_i|<^;_i), 6^ = 0_, and = 0. Ob- 

viously, the representation matrix [T-^]/.;/ = ((^/|H^|(y9;/) of is tridiagonal in the 
L-dimensional Hilbert space spanned by the {\'Pi)}i=o,...,L-i {L ^ D): 

... \ 

... 

... 

64 ... 

■■. / 

Applying the Lanczos recursion (|ll|l . the eigenvalues En and eigenvectors \^n) of 
are approximated by 



(12) 



/ ao 61 

bi ai 62 

62 0,2 bs 

63 03 

V ; ; ; ; 



(13) 



L-l 

En and \i^t) = Y.cii\^i), 

1=0 
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respectively, where the L coefficients ; are the components of the {n — th) eigenvectors 
of with eigenvalue Ef^ . The eigenvalue spectrum of can be easily determined using 
standard routines from libraries such as EISPACK (see http:/ /www. nctlib.org). Increas- 
ing L we check for the convergence of an eigenvalue of in a specific energy range. So 
we can avoid spurious eigenvalues for fixed Lanczos dimension L which disappear as one 
varies L 

Note that the convergence of the Lanczos algorithm is excellent at the edges of the 
spectrum (the ground state for example is obtained with high precession using at most 
^ 100 Lanczos iterations) but rapidly worsens inside the spectrum. So Lanczos is suitably 
used only to obtain the ground state and a few low lying excited states. 

2'3.2. Implementation of matrix vector multiplication . The core operation of most 
ED algorithms is a MVM. It is quite obvious that our matrices are extremely sparse 
because the number of non-zero entries per row of our Hamilton matrix scales linearly 
with the number of electrons. Therefore a standard implementation of the MVM step 
uses a sparse storage format for the matrix, holding the non-zero elements only. Two 
data schemes are in wide use, the compressed row storage (CRS) and the jagged diagonal 
storage (JDS) format ^31, where the latter is the method of choice for vector computers. 
The typical storage requirement per non-zero entry is 12-16 Byte for both methods, i.e. 
for a matrix dimension of Z? = 10^ about one TByte main memory is required to store 
only the matrix elements of the EP Hamiltonian. Both variants can be applied to any 
sparse matrix structure and the MVM step can be done in parallel by using a parallel 
library such as PETSc (see http://www-unix.mcs.anl.gov/petsc/petsc-as/). 

To extend our EP studies to even larger matrix sizes we store no longer the non-zero 
matrix elements but generate them in each MVM step. Of course, at that point standard 
libraries are no longer useful and a parallel code tailored to each specific class of Hamil- 
tonians must be developed. For the Holstein-Hubbard EP model we have established a 
massively parallel program using the Message Passing Interface (MPI) standard. The 
minimal total memory requirement of this implementation is three vectors with Hilbcrt 
space dimension. 

The parallelization approach follows the inherent natural parallelism of the Hilbert 
space, which can be constructed as the tensorial product space of electrons and phonons 
{\b) = |e) (g) |p)} (cf. subsubscct . 12 ' 1 . l|l . Assuming, that the electronic dimension (Dg) is 
a multiple of the number of processors used (A'cpu) we can easily distribute the electronic 
basis states among these processors, i.e. processor i{0 < i < iVcpu— 1) is holding the basis 
states (ci = iDe/Ncpu -I- 1, . . . , (i -I- l)£'e/7Vcpu)- As a consequence of this choice only the 
electronic hopping term generates inter-processor communication in the MVM while all 
other (diagonal electronic) contributions can be computed locally on each processor. 

Furthermore, the communication pattern remains constant within a single run for 
all MVM steps and the message sizes (at least Dp words) are large enough to ignore 
the latency problems of modern interconnects. Using supercomputers with hundreds of 
processors and one TBytes of main memory, such as IBM p690 clusters or SGI Altix 
systems, we are able to run simulations up to a matrix dimension of 30 x 10^. 
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2'4. Algorithms for estimating spectral functions. - The numerical caleulation of spec- 
tral functions 



A'^iu}) = - lim -Im 
b-1 



(V-olOt O|^o) 

uj — ti. + + 11] 



(14) = \{i^nm^)\H[u-{E,,-Eo)l 

?i=0 

where O is the matrix representation of a certain operator O (e.g., the creation operator 
of an electron with wave number k if one wants to calculate the single-particle spectral 
function; or the current operator j = ~ '^I+i'^i) '^'^'^ interested in the 

optical conductivity), involves the resolvent of the Hamilton matrix H. Once we have 
obtained the eigenvalues and eigenvectors of H we can plug them into eq. I|14|) and 
obtain directly the corresponding dynamical correlation or Green functions. In practice 
this 'naive' approach is applicable for small Hilbert spaces only, where the complete 
diagonalization of the Hamilton matrix is feasible. 

For the typical EP problems under investigation we deal with Hilbert spaces having 
total dimensions D of 10^-10^^. Finding all eigenvectors and eigenstates of such huge 
Hamilton matrices is impossible, because the CPU time required for exact diagonalization 
of H scales as and memory as . Fortunately, there exist very accurate and well- 
conditioned linear scaling algorithms for a direct approximate calculation of A'-'{uj). 

2'4.1. Lanczos recursion method. Having determined the ground state ItAo ) by the 
Lanczos technique, we can use again the recursion relation but with the initial 

state |(/3o) = OjV'Q )/ ■\/ (t/'q jOl'OI'^Q ), to determine within the so-called Lanczos recursion 
method (LRM) or spectral decoding method (SDM) an approximative spectral function, 

L-l 

(15) A^(c.) = KoP(V'olOtOIV'o) 5[u^ ~ {E^ - E^)l 
or equivalently 

(16) A-(.) = - lim ilm (^olOtOl^o) 



uj + i?7 — ao 

w -\-\ri — ai — ■ 



?)2 
'^2 



02 

which is built up by L 5-peaks. 

Of course, the true spectral function A'-' {lo) has D (5-peaks. According to the Lanczos 
phenomenon, the approximated spectral weights and positions of the peaks converge 
to their true values with increasing L. Some of the main problems of the LRM/SDM 
are: (i) The convergence is not uniform in the whole energy range, (ii) There exist 
so-called spurious peaks, which appear and disappear as L is increased, i.e., when the 
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iteration proceeds, (iii) Without computationally expensive re-orthogonalization only a 
few hundred iterations are possible. 

2'4.2. Kernel polynomial method. The idea behind a conceptionally different ap- 
proach, the kernel polynomial method (KPM) (for a review see jl4|'). is to expand A'^{lo) 
in a finite series of L + 1 Chebyshev polynomials Tm{x) = cos[marccos(x)]. Since the 
Chebyshev polynomials are defined on the real interval [—1,1], we apply first a sim- 
ple linear transformation to the Hamiltonian and all energy scales: X = (H — b)/a, 
X = {lu — b)/a, a = (i?max ^ £'inin)/2(l — e), and b = (E'inax + ^'min)/2 (the small constant 
e is introduced in order to avoid convergence problems at the endpoints of the interval 
-a typical choice is e ^ 0.01 which has only 1% impact on the energy resolution [ISp. 
Then the expansion reads 

(17) A°{x) = -jL= + 2 X] t^ZTmix)] , 
with the coefficients (moments) 

(18) " / 1 ^^^'"(^)^°(^) = (V^o|OtT™(X)0|V'o). 
Equation H17|) converges to the correct function for L ^ oo. Again the moments 

(19) ^J'2m = '2{(f>.m\(l),n) - and H2m+1 = 2((/),„+i |0,„) - ^if 

can be efficiently obtained by repeated parallelized MVM, where \(j),n+i) = 2X|(/)m) — 
|0m-i) but now \(f>i) — X|0o) and \(f>o) — 0\tJjo) with j^'o) determined by Lanczos ED. 

As is well known from Fourier expansion, the series (|17|l with L finite suffers from 
rapid oscillations (Gibbs phenomenon) leading to a poor approximation to A'~'{uj). To 
improve the approximation the moments //„ are modified fin dnfJ-m where the damping 
factors gn are chosen to give the 'best' approximation for a given L. This modification is 
equivalent to a convolution of the infinite scries with a smooth approximation K^ix^y) 
to 5{x — y), a so-called approximation kernel. The appropriate choice of this kernel, that 
is of gn, e.g. to guarantee positivity of A^{uj), lies at the heart of KPM. We mainly use 
the Jackson kernel which results in a uniform approximation whose resolution increases 
as 1/i, but for the determination of the single-particle Green functions below we use a 
Lorentz kernel which mimics a finite imaginary part rj in eq. H14|l . see |14| . 

In view of the uniform convergence of the expansion, KPM is a method tailored to 
the calculation of spectral properties. Most important, spectral functions obtained via 
KPM are not subject to uncontrolled or biased approximations: The accuracy of its 
outcome depends only on the expansion depth L, and can be made as good as required 
by just increasing L. Of course one is restricted to finite systems of moderate size whose 
associated Hamilton matrix does not exceed available computational resources. 
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2'4.3. Cluster perturbation theory (CPT). The spectrum of a finite system of sites 
which we obtain through KPM differs in many respects from that in the thermodynamic 
hmit N ^ oo, cspeciahy it is obtained for a finite number of momenta K = nm/N only. 
The most obvious feature is the resulting discreteness of energy levels which is a property 
already of the non-interacting system. While we cannot easily increase N without going 
beyond computationally accessible Hilbert spaces, we can try to extrapolate from a finite 
to the infinite system. 

For this purpose we first calculate the Green function (oj) for all sites i,j = 1, . . . ,N 
of a A'^-sizc cluster with open boundary conditions, and then recover the infinite lattice 
by pasting identical copies of this cluster at their edges. The 'glue' is the hopping V 
between these clusters, where Vki = t for \k — l\ ~ 1 and k,l = 0,1 mod N, which is dealt 
with in first order perturbation theory. Then the Green function Gtj {to) of the infinite 
lattice is given through a Dyson equation 

(20) G^.iLu) = Gt.H + J2 Gtkiio)VkiGi,iLu), 



where indices of G'^(tt') are counted modulo N . Obviously this order of perturbation in 
V is exact for the non-interacting system. We thus get rid of the discreteness addressed 
above. The Dyson equation is solved by Fourier transformation over momenta K = kN 
corresponding to translations by N sites 



(21) G,,iK,iu) 



G^iiu) 



1 - V{K)G''{uj) 



from which one finally obtains 



(22) G{k,co) - ^ E G'i,{Nk,Lo)cM-^H^-j))- 

In this way, which is called CPT |16| . we obtain a Green function G{k, uj) with contin- 
uous momenta k from the Green function G'lj {uj) on a finite cluster. Two approximations 
are made, one by using first order perturbation theory in V = t, the second on assuming 
translational symmetry in Gij (lu) which is only approximatively satisfied. 

In principle, the CPT spectral function G{k, uj) does not contain any more information 
than the cluster Green function G^j{u>) already does. But extrapolating to the infinite 
system it gives a first hint at the scenario in the thermodynamic limit. However, CPT 
does not describe effects which only occur on large length scales, like Anderson local- 
ization (see the paper by Fehske, Bronold and Alvermann) or the critical behavior at a 
phase transition. Providing direct access to spectral functions, still without relying on 
possibly erroneous approximations, CPT occupies a niche between variational approaches 
like (D)DMRG (see sect. 151 and 1^ and methods directly working in the thermodynamic 
limit like the variational ED method j9). 
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3. Density matrix renormalization group approach 

The Density Matrix Renormalization Group (DMRG) is one of the most powerful nu- 
merical techniques for studying many-body systems. It was developed by Steve White |17| 
in 1992 to overcome the problems arising in the application of the standard Numerical 
Renormalization Group (NRG) to quantum lattice many-body systems such as the Hub- 
bard model Since then the approach has been extended to a great variety of problems, 
from Classical Statistical Physics to Quantum Chemistry (including ab initio calculations 
of electronic structures in molecules) and, recently, to Nuclear Physics and the Physics of 
Elementary Particles and Fields. A review article on DMRG and its numerous applica- 
tions has recently been published |18| . Additional information can also be found on the 
DMRG web page at http://www.dmrg.info. A detailed discussion of the basic DMRG al- 
gorithms and their implementation has been published in ref. |19| . Readers interested in 
a simple example should consider the application of DMRG to single-particle problems, 
which is also discussed there. The source code of a single-particle DMRG program (in 
the programming language C'''^) is now part of the ALPS distribution |2U| . 

DMRG techniques for strongly correlated systems have been substantially improved 
and extended since their conception and have proved to be both extremely accurate 
for low-dimensional problems and widely applicable. They enable numerically exact 
calculations (as good as exact diagonalizations) of low-energy properties on large lattices 
with up to a few thousand particles and sites (compared to less than a few tens for exact 
diagonalizations). The calculation of high-energy excitations and dynamical spectra for 
large systems has proved to be more difficult and will be discussed in sect. HI 

3'1. Renormalization group and density matrix. - Consider a quantum lattice system 
with TV sites (in general, we are interested in the case ^ oo or at least N ^ 1). 
The Hilbert space of this system is the Fock space of all properly symmetrized many- 
particle wave functions. As seen in sect. |51 its dimension D grows exponentially with 
the number of sites N. Obviously, a lot of these states are not necessary to investigate 
specific properties of a model such as the ground state. Thus, a number of methods have 
been developed to perform a projection onto a subspace of dimension d <^ D and then 
an exact diagonalization of the Hamiltonian in this subspace. 

Such an approach, called the Numerical Renormalization Group (NRG), was devel- 
oped by Wilson 30 years ago to solve the Kondo impurity problem f2]- The key idea 
is the decomposition of the system into subsystems with increasing size (see the paper 
on the NRG method by Hewson). The subsystem size is increased by one site at each 
step as shown in fig.^ Each subsystem is diagonalized successively and the information 
obtained is used to truncate the Hilbert space before proceeding to the next larger sub- 
system. Let m and n be the dimension of the (effective) Hilbert spaces associated with 
the subsystem made of the first £ sites and with the site £ + 1, respectively. A basis of 
dimension d = mn for the next subsystem with £ + I sites is built as a tensor product of 
the subsystem and site bases. Assuming that d is small enough, the effective Hamiltonian 
can be fully diagonalized in the tensor-product basis. The energy is used as a criterion 
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Fig. 4. - Schematic representations of the usual NRG algorithm (left), an alternative NRG 
algorithm (middle) and the energy levels in the effective Hamiltonian He for increasing subsystem 
system size £ (right). 



to truncate the subsystem Hilbert space. The lowest m eigenstates are kept to form 
a new basis while the high-energy eigenstates are discarded. The new subsystem with 
£ + I sites and an effective Hilbert space of dimension m can be used to start the next 
iteration. In summary, this procedure provides a transformation = R[Hi\ forming 

effective Hamiltonians of fixed dimension m which describe the low-energy physics of 
increasingly larger systems. In this transformation the high-energy states are steadily 
traced out as the system grows as illustrated in fig. 01 Such a transformation R is usually 
called a renormalization group (RG) transformation. A different implementation of the 
NRG idea is possible if the system is homogeneous like in the Hubbard model. A copy 
of the current subsystem can be substituted for the added site in the above procedure. 
Thus the subsystem size doubles at every iteration as illustrated in fig. 0] 

The NRG method is very accurate for the Kondo problem and more generally for 
quantum impurity problems. Unfortunately, NRG and related truncation schemes have 
proved to be unreliable for quantum lattice systems such as the Hubbard model \'22\ . 
It is easy to understand the failure of the standard NRG in those cases. A subsystem 
always has an artificial boundary at which the low-energy eigenstates of a quantum 
lattice Hamiltonian tend to vanish smoothly. Thus the truncation procedure based on 
the effective cigenencrgies may keep only states that vanish at the artificial boundary. As 
a consequence, at later RG iterations the eigenstates of the effective Hamiltonian of larger 
subsystems may have unwanted features like nodes where the artificial boundary of the 
previous subsystems were located. The application of the second NRG algorithm (in the 
middle in fig. ^ to the problem of a quantum particle in a one-dimensional box gives a 
clear illustration of this effect |19H23| . The ground state wavefunction for a TV-site lattice 

^i^) ~ \J n\i ^^^^ ( jv+i ) ^^^^ ^ minimum (node) where the ground state wavefunction 
of the twice larger system {N 2N) has a maximum as seen in fig. Therefore, 
the low-energy eigenstates of a small system are not necessarily the best states to form 
the low-energy eigenstates of a larger system. An approach proposed by White and 
Noack [52! to solve this problem is the construction of an effective Hamiltonian including 
the effects of the subsystem environment to eliminate the artificial boundary. DMRG is 
the extension of this idea to interacting many-particle problems. 
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Fig. 5. - Ground-state wavefunctions of the tight-binding particle-in-the-box problem for two 
systems of A'' =10 sites (circles) and one system of A'^ = 20 sites (squares). 



Consider a quantum system which can be spht into two parts called the subsystem 
and its environment. A basis of the system Hilbcrt space is given by the tensor product 

(23) N,j) = N)®|j)e 

of basis states \i) (i = l,...,d) and |j)E (j = Ij-.-jcJe) for the subsystem and its 
environment, respectively. The dimension of this basis is ds = d ■ d-E- Any state \tp) of 
the system can be expanded in this basis \iJj) = Si j '/'u N) b)E • The most important 
states in the subsystem to represent the state {ip) are given by its reduced density matrix, 
which is obtained by tracing out the states of the environment 

(24) p.,.. 

This density matrix is symmetric and has d positive eigenvalues Wa > satisfying the 
normalization '^^Wa = 1. A new basis of the subsystem, \va) = '^^Vai\i), can be 
defined using the eigenvectors of the density matrix pi.i'Vai' — WaVai] a,i = 1, . . . , d. 
In the new basis the state \Tp) can be written 

(25) W ^^Xa\Va) ®\Ua)E, 

a 

with = Wa > and normalized states |uq)e = j ''^ai''Pij\j)'s^- Therefore, Wa is 

the probability that a subsystem is in a state \va) if the superblock is in the state j?/;). 
The density matrix provides an optimal choice for selecting the m states of the subsystem 
to be kept in a RG transformation: keep density matrix eigenstates \va) with the largest 
weights Wa ■ This is the key idea of the density matrix renormalization group approach. 

3'2. DMRG algorithms . - In a DMRG calculation one first forms a new subsystem 
with i+1 sites and an effective Hilbert space of dimension d = mn by adding a site to the 
current subsystem with (. sites as in the NRG method. Then one considers a larger system, 
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called a superblock, which is made of the new subsystem and an environment. A basis 
of the superblock Hilbert space is given by the tensor product (|23|l . Assuming initially 
that wc want to compute the system ground state only, wc then calculate the ground 
state of the superblock Hamiltonian using the techniques discussed in subsect. 12 31 
Then the reduced density matrix H24|l for the subsystem is obtained by tracing out the 
environment states. The m density matrix eigenstates \va) with the largest weights Wa 
are kept to build an optimal effective basis of dimension m for the subsystem with £ +1 
sites. This subsystem is then used to start the next iteration and build the next larger 
subsystem. This procedure defines a generic density matrix RG transformation. Clearly, 
the accuracy of a DMRG calculation will depend on the quality of the environment used. 
The environment should make up as much as possible of the lattice which is not already 
included in the subsystem. Constructing such a large and accurate environment is as 
difficult as the original quantum many-body problem. Therefore, the environment must 
also be constructed self-consistently using a density matrix RG transformation. 

In his initial papers jJTj, White described two DMRG algorithms: the infinite-system 
method and the finite-system method. The infinite-system method is certainly the sim- 
plest DMRG algorithm and is the starting point of many other algorithms. Its defining 
characteristic is that the environment is constructed using a "reflection" of the current 
subsystem. The superblock size increases by two sites at each step as illustrated in fig. 
The system is assumed to be homogeneous and "symmetric" to allow this operation. It- 
erations are continued until an accurate approximation of an infinite system is obtained. 
The infinite-system method is sometimes used to calculate properties of finite systems. 
While this approach may work just fine in many cases, it should be kept in mind that 
there is no guarantee of convergence to the eigenstates of the finite system. 

The finite-system method is the most versatile and reliable DMRG algorithm and with 
some later improvements it has also become a very efficient method. It is designed 

to calculate the properties of a finite system accurately. The environment is chosen so 
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Fig. 6. - Schematic representation of DMRG algorithms. Left: Infinite-system algorithm (from 
top to bottom). Right: Finite-system DMRG algorithm for a ten-site lattice. The environment 
block is on the right-hand side when going from top to bottom and on the left-hand-side when 
going from bottom to top. 
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that the superblock represents the full lattice at every iteration. The environments can 
also be considered as being subsystems and are calculated self-consistently using DMRG 
with the usual subsystems playing the part of their environment. Iterations arc continued 
back and forth through every configuration of the superblock (i.e., every splitting of the 
lattice in one subsystem and its environment) until convergence. This procedure is 
illustrated in fig. This ensures the self-consistent optimization of the subsystems and 
their environments (for a given number m of states kept) and thus considerably improves 
the results reliability compared to the infinite-system method. 

Most DMRG algorithms use only two blocks (the subsystem and its environment) 
to form the superblock, but it is possible and useful for some problems to consider a 
more complicated configuration. For instance, one can use four blocks to treat one- 
dimensional problems with periodic boundary conditions and using several blocks can 
also be advantageous for systems with boson degrees of freedom such as phonons |24j . 

The DMRG sites usually correspond to the physical sites of the lattice model in- 
vestigated, such as spin, fermionic, or bosonic degree of freedom. [For a phonon mode 
(boson), which has an infinite Hilbcrt space, a DMRG site represents a finite dimensional 
basis of the phonon states as already discussed for exact diagonalization methods 
fsubsect. I2TII .] However, the DMRG sites can also represent a combination of several 
physical sites [for instance, the electron and phonon at each site of the Holstein-Hubbard 
model (0)], a fraction of the Hilbert space associated with a given physical site (as in the 
pseudo-site method for bosonic degrees of freedom presented in subsubsect. I3"4.l1l . 

It is possible to compute several quantum states simultaneously with DMRG. In that 
case, the density matrix is formed as the sum of the density matrices (|24|l calculated for 
each target state. A target state can be any quantum state which is well-defined (and 
can be computed) in every superblock of a DMRG calculation. This feature turns out to 
be very important for the calculations of dynamical properties (see sect. 

3'3. Truncation errors. - With DMRG an error is made when projecting operators 
onto the subspace spanned by the most important m density matrix eigenstates. This 
is called the truncation error. Probably the most important characteristic of a DMRG 
calculation is the rate at which the truncation error decreases with an increasing number 
m of states kept. In the most favorable cases (gapped one-dimensional systems with 
short-range interactions only and open boundary conditions), the accuracy increases 
roughly exponentially with m. For instance, the ground-state energy of the spin-one 
Heisenberg chain on lattices with hundreds of sites can be calculated to an accuracy of the 
order of 10~^° with a modest computational effort. In very difficult cases (long-range off- 
diagonal interactions, two-dimensional systems with periodic boundary conditions), the 
truncation error in the ground state energy can decrease as slowly as m~^. It is possible 
to calculate exactly the density matrix spectrum of several integrable models |26| . This 
analysis shows that the distribution of density matrix eigenvalues Wa varies greatly. As 
a result, truncation errors may fall exponentially with increasing m in some favorable 
cases but decrease extremely slowly for other ones. Correspondingly, DMRG accuracy 
and performance depend substantially on the specific problem investigated. 
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For any target state \Tp) written down in its representation H25|l . the truncation of the 
density matrix basis corresponds to making an approximation 

(26) = ^ A„ ® |u„)e, 

m kept states 

which minimizes the error 

(27) ^ \[^) - \^/)f = Wa^l- 

d—rn discarded states m kept states 

for a fixed number m of states kept. The total weight of the discarded density matrix 
eigenstates (discarded weight) Dm is related to the truncation errors of physical quan- 
tities. For Dm ^ 1 it can be shown that the truncation error in the eigenenergy of 
any target eigenstates of the Hamiltonian H scales linearly with the discarded weight, 
E'm - Sexact = cDm + 0(i?^) , whcrc = |77 1^') IS the energy in the 
approximate state and c is a constant. For the expectation values of other operators 
the truncation error theoretically scales as \/Dm- The discarded weight Dm decreases 
(and thus the accuracy of DMRG results increases) when the number m of density ma- 
trix eigenstates kept is increased. In particular, for large enough m the discarded weights 
vanish at every RG iteration for both subsystems and environments and truncation errors 
become negligible. Therefore, DMRG is an exact numerical method as defined in the 
introduction (sect.^. Moreover, if the number m of density matrix eigenstates kept is 
so large that the discarded weight is exactly zero at every RG iteration, DMRG becomes 
equivalent to an exact diagonalization (scct.[5J). 

In real DMRG applications, series of density matrix basis truncations are performed 
in successive supcrblock configurations. Therefore, the measured discarded weights Dm 
do not represent the real error in the wave function of the target. Nevertheless, in most 
cases, the energy truncation error scales linearly with the average measured discarded 
weight Dm for Dm ^ 1 as shown in fig. IHa). This can (and should) be used to make 
an Dm ^ extrapolation of the energy. For other physical quantities such as an order 
parameter, truncation errors sometimes scale as {DmY , with r « 0.5 as seen in fig.Cfb). 
This can also be used to estimate truncation errors. The measured discarded weight Dm 
alone should not be used as an indication of a calculation accuracy, because truncation 
errors for physical quantities can be several orders of magnitude larger than Dm- 

It should be noted that DMRG can be considered as a variational approach. The 
system energy E^i/j) = {ip\H\ip) / {ip\ip) is minimized in a variational subspace of the 
system Hilbert space to find the ground-state wavefunction jV'o) smd energy Eq = E{tpQ). 
If the ground-state wavefunction is calculated with an error of the order of e ~ \/Dm ^ 1 
(i.e., IV') = iV'o) + £10); with {(\)\(\)) = 1), the energy obtained is an upper bound to the 
exact result and the error in the energy is of the order of (~ Dm) as in all variational 
approaches. For specific algorithms the variational wave function can be written down 
explicitly as a matrix product wave function or as a product of local tensors |18| . 
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Fig. 7. - (a) Ground state energy E'^ calculated with DMRG for a 12 x 3 Hubbard ladder with 
U = t and 6 holes as a function of the discarded weight Dm for several number of density matrix 
eigenstates kept m = 300 to 2200. The zero of the energy is given by a (nearly) exact quantum 
Monte Carlo result 1271 . (b) Staggered bond order parameter 5m calculated with DMRG in an 
extended ID half-filled Hubbard model I28| for U = 3t and V = 1.5t on a 1024-site lattice as a 
function of V Dm for m = 200 to 1200. Solid lines are linear fits. 



The computational effort of a DMRG calculation usually increases as a power law 
for increasing system size N, number m of states kept, or DMRG site dimension n. In 
the most favorable case (ID system with short-range interactions only), the CPU time 
increases theoretically as Nm^n^, while the amount of stored data is of the order of 
~ m^n^. In most cases, however, m has to be increased with the system size TV to keep 
truncation errors constant. 

3'4. Methods for electron-phonon systems. - A significant limitation of DMRG and 
exact diagonalizations is that they require a finite basis for each site. In electron-phonon 
lattice models such as the Holstein-Hubbard model (0) , the number of phonons (bosons) 
is not conserved and the Hilbert space is infinite for each site representing an oscillator. 
Of course, the number of phonons can be artificially constrained to a finite number M 
per site, but the number M needed for an accurate treatment may be quite large. This 
often severely constrains the system size or the regime of coupling which may be studied 
with DMRG [^S] and exact diagonalizations (sect. l2)l. 

Here, we describe two methods for treating systems including sites with a large Hilbert 
space. Both methods use the information contained in a density matrix to reduce the 
computational effort required for the study of such systems. The first method called 
pseudo-site method, is just a modification of the usual DMRG technique which allows us 
to deal more efficiently with sites having a large Hilbert space. The second method |3(JI 
I31| . called the optimal basis method, is a procedure for generating a controlled truncation 
of the phonon Hilbert space, which allows the use of a very small optimal basis without 
significant loss of accuracy. 
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Fig. 8. - Symbolic representations (a) of the standard DMRG approach for Af = 8 and (b) of 
the pseudo-site approach for P = 3. 



3'4.1. Pseudo-site method. The DMRG algorithms presented above can easily be 
generalized to treat systems including phonons (bosons). In a standard implementation 
of the DMRG method, however, each boson forms one lattice site and thus memory and 
CPU time requirements increase as and M^, respectively. Therefore, performing 
calculations for the Holstein model requires much more computer resources than compu- 
tations for purely electronic systems. 

To understand the basis of the pseudo-site approach j29) . it is important to note 
that, in principle, the computer resources required by DMRG increase linearly with the 
number of lattice sites. Thus, DMRG is much better able to handle several few-state 
sites rather than one many-state site. The key idea of the pseudo-site approach is to 
transform each boson site with M — 2^ states into P pseudo-sites with 2 states. This 
approach is motivated by a familiar concept: The representation of a number in binary 
form. In this case the number is the boson state index s going from to P — 1. Each 
binary digit rj is represented by a pseudo-site, which can be occupied [vj = 1) or empty 
{rj = 0). One can think of these pseudo-sites as hard-core bosons. Thus, the level (boson 
state) with index s = is represented by P empty pseudo-sites, while the highest level, 
s = 2^ — 1, is represented by one boson on each of the P pseudo-sites. 

Figure |H1 illustrates the differences between standard and pseudo-site DMRG ap- 
proaches for M = 8 (P = 3). In the standard approach [fig. |HIa)], a new block (dashed 
rectangle) is built up by adding a boson site (oval) with M states to another block (solid 
rectangle) with m states. Initially, the Hilbert space of the new block contains niM 
states and is truncated to m states according to the DMRG method. In the pseudo-site 
approach [fig. Elb)], a new block is made of the previous block with m states and one 
pseudo-site with two states. The Hilbert space of this new block contains only 2m states 
and is also truncated to m states according to the DMRG method. It takes P steps to 
make the final block (largest dashed rectangle) including the initial block and all pseudo- 
sites, which is equivalent to the new block in fig. |HIa). However, at each step we have to 
manipulate only a fraction 2/M of the bosonic Hilbert space. 

To implement this pseudo-site method, we introduce P pseudo-sites j ~ 1, P with a 
two-dimensional Hilbert space {[j'j), fj = 0, 1} and the operators aj, such that aj\l) = 
|0), fljlO) = and is the hermitian conjugate of a^. These pseudo-site operators 
have the same properties as hard-core boson operators: cljOj + ajUj = 1, and operators 
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on different pseudo-sites commute. The one-to-one mapping between a boson level |s), 
s ~ 0, M —1, where b^b\s) = s\s), and the P-pseudo-site state |ri, r2, rp) is given by 
the relation s — X^^Li between an integer number and its binary representation. 

The next step is to write all boson operators in terms of pseudo-site operators. It 
is obvious that the boson number operator is given by iVb = b^b = X]jLi ^-'"^ aja^. 
Other boson operators take a more complicated form. For instance, to calculate the 
representation of b^ we first wr ite b^ = Bt/ZYbTT, where B^\s) = \s + l). The pseudo- 
site operator representation of the second term is 

(28) V^^b + l = J2 ^i(^i) A2{r2)...Ap{rp), 

s=Q 

where Aj{l) = a]a^, Aj{0) = a^a] and the Vj {j = 1, .., P) are given by the binary digits 
of s. For Pt 

we find 

(29) P^ = a| + fljOj^ + 03030]^ + ... + apflp-^apj-'-Q-i- 

Thus one can substitute P = log2 (M) pseudo-sites for each boson site in the lattice and 
rewrite the system Hamiltonian and other operators in terms of the pseudo-site operators. 
Then the finite system DMRG algorithm can be used to calculate the properties of this 
system of interacting electrons and hard-core bosons. 

The pseudo-site approach outperforms the standard approach when computations be- 
come challenging. For M = 32, the pseudo-site approach is already faster than the stan- 
dard approach by two orders of magnitude. With the pseudo-site method it is possible to 
carry out calculations on lattices large enough to eliminate finite size effects while keep- 
ing enough states per phonon mode to render the phonon Hilbert space truncation errors 
negligible. For instance, this technique provides some of the most accurate results |29II32| 
for the polaron problem in the Holstein model (see also the paper by Fehske, Alvermann, 
Hohcnadler and Wellcin). It has also been successfully used to study quantum phase 
transitions in the ID half-filled Holstcin-Hubbard model (see ref. [33] and the separate 
paper by Fehske and Jeckelmann). 

3'4.2. Optimal phonon basis. The number of phonon levels M needed for an accurate 
treatment of a phonon mode can be strongly reduced by choosing a basis which minimizes 
the error due to the truncation of the phonon Hilbert space instead of using the bare 
phonon basis made of the lowest eigenstates of the operators b\b^. As with DMRG, in 
order to eliminate phonon states without loss of accuracy, one should transform to the 
basis of eigenvectors of the reduced density matrix and discard states with low probability. 
The key difference is that here the subsystem is a single site. To be specific, consider the 
translationally invariant Holstein-Hubbard model A site includes both the phonon 
levels and the electron degrees of freedom. Let a label the four possible electronic states 
of a particular site and let s label the phonon levels of this site. Let j label the combined 
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states of all of the rest of the sites. Then a wave function of the system can be written 
as 



where r also labels the phonon levels of this site. Let Wak be the eigenvalues and 4>ak{n) 
the eigenvectors of p, where k labels the different eigenstates for a given electronic state of 
the site. The Wak are the probabilities of the states 4>ak if the system is in the state (|30ll . 
If Wak is negligible, the corresponding eigenvector can be discarded from the basis for the 
site without affecting the state (|3UII . If one wishes to keep a limited number of states m 
for a site, the best states to keep are the eigenstates of the density matrices with the 
largest eigenvalues. In EP systems, these m eigenstates form an optimal phonon basis. 
We note that we obtain different optimal phonon states for each of the four electron 
states of the site. 

Unfortunately, in order to obtain the optimal phonon states, we need the target 
state l|5n)l. which we do not know -usually we want the optimal states to help get this 
state. This problem can be circumvented in several ways PU]. Here, we describe one 
algorithm in conjunction with an exact diagonalization approach (sect.|51) but it can also 
be incorporated into a standard DMRG algorithm. One site of the system (called the 
big site) has both optimal states and a few extra phonon levels. (These n extra levels are 
taken from a set of M ^ m bare levels but are explicitly orthogonalized to the current 
optimal states.) To be able to perform an exact diagonalization of the system, each site 
of the lattice is allowed to have only a small number of optimal phonon levels, m ~ 3 — 4. 
This approach is illustrated in fig. |21 The ground state of the Hamiltonian is calculated 
in this reduced Hilbert space using an exact diagonalization technique. Then the density 
matrix 131() of the big site is diagonalized. The most probable m eigenstates are new 
optimal phonon states, which are used on all other sites for the next diagonalization. 
Diagonalizations must be repeated until the optimal states have converged. Each time 
different extra phonon levels are used for the big site. They allow improvements of the 
optimal states by mixing in the M bare states little by little. 

The improvement coming from using optimal phonon states instead of bare phonon 
levels is remarkable. For the Holstein-Hubbard model the ground state energy converges 
very rapidly as a function of the number of optimal phonon levels |3UI I31| . Two or three 
optimal phonon states per site can give results as accurate as with a hundred or more 
bare phonon states per site. For intermediate coupling (wq = f , g = 1.5, and U ~ Q) in 
the half-filled band case, the energy is accurate to less than 0.1% using only two optimal 
levels, whereas keeping eleven bare levels the error is still greater than 5%. 



(30) 




The density matrix for this site for a given electronic state a of the site is 



(31) 
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Fig. 9. - Schematic representation of the big site algorithm. Each site (circles) has the electronic 
degrees of freedom and three optimal states (wiggly bars). The big site (second site from the 
left) has the optimal states plus two bare levels (straight bars). 



Combined with ED techniques the above algorithm for generating optimal phonon 
states can be significantly improved |34[I35) . First, the use of a different phonon basis for 
the big site artificially breaks the system symmetries. This can be solved by including all 
those states into the phonon basis that can be created by symmetry operations and by 
summing the density matrices generated with respect to every site. Second, the effective 
phonon Hilbert space is unnecessarily large if one uses the configuration shown in fig. |51 
As eigenvalues Wak of the density matrix decrease very rapidly we can introduce a cut- 
off for the lattice phonon states, which is reminiscent of the energy or phonon number 
cut-off (jHl discussed in subsubscct. 12 1.21 to further reduce the phonon Hilbert space 
dimension without loss of accuracy. 

In the Holstein-Hubbard model it is possible to transfer optimal phonon states from 
small systems to larger ones because of the localized nature of the phonon modes. There- 
fore, one can first calculate a large optimal phonon basis in a two-site system and then 
use it instead of the bare phonon basis to start calculations on larger lattices. This is the 
simplest approach for combining the optimal phonon basis method with standard DMRG 
techniques for large lattices such as the infinite- and finite-system algorithms pj>6] . 

The features of the optimal phonon states can sometimes be understood qualita- 
tively ISOl 1^ . In the weak-coupling regime optimal states are simply eigenstates of an 
oscillator with an equilibrium position (g) Ri 2g as predicted by a mean-field approxima- 
tion. In the strong-coupling regime {g^oJo ^ U, t) the most important optimal phonon 
states for = or 2 electrons on the site can be obtained by a unitary Lang-Firsov 
transformation of the bare phonon ground state S{g) = e"^^*^^'"''*-'"' , in agreement 
with the strong-coupling theory. The optimal phonon state for a singly occupied site 
is not given by the Lang-Firsov transformation but is approximately the superposition 
of the optimal phonon states for empty or doubly occupied sites due to retardation ef- 
fects |3U| . In principle, the understanding gained from an analysis of the optimal phonon 
states calculated numerically with our method could be used to improve the optimized 
basis used in variational approaches (see the related paper by Cataudella). 

An interesting feature of the optimal basis approach is that it provides a natural 
way to dress electrons locally with phonons \6\\ . This allows us to define creation and 
annihilation operators for composite electron-phonon objects like small polarons and 
bipolarons and thus to calculate their spectral functions. However, the dressing by 
phonons at a finite distance from the electrons is completely neglected with this method. 
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4. Dynamical DMRG 

Calculating the dynamics of quantum many-body systems such as solids has been 
a long-standing problem of theoretical physics because many experimental techniques 
probe the low-energy, low-temperature dynamical properties of these systems. For in- 
stance, spectroscopy experiments, such as optical absorption, photoemission, or nuclear 
magnetic resonance, can measure dynamical correlations between an external perturba- 
tion and the response of electrons and phonons in solids |37| . 

The DMRG method has proved to be extremely accurate for calculating the proper- 
ties of very large low-dimensional correlated systems and even allows us to investigate 
static properties in the thermodynamic limit The calculation of high-energy excitations 
and dynamical spectra for large systems has proved to be more difficult and has be- 
come possible only recently with the development of the dynamical DMRG (DDMRG) 
approach Here we first discuss the difficulty in calculating excited states and 

dynamical properties within the DMRG approach and several techniques which have 
been developed for this purpose. Then we present the DDMRG method and a finite-size 
scaling analysis for dynamical spectra. In the final section, we discuss the application of 
DDMRG to electron-phonon systems. 

4'1. Calculation of excited states and dynamical properties. - The simplest method 
for computing excited states within DMRG is the inclusion of the lowest R eigenstates 
as target instead of the sole ground state, so that the RG transformation produces 
effective Hamiltonians describing these states accurately. As an example, fig. ^| shows 
the dispersion of the lowest 32 eigenenergies as a function of the momentum Hk in the 
one-dimensional Holstein model on a 32-site ring with one electron (the so-called polaron 
problem, see the paper by Fehske, Alvermann, Hohenadler, and Wellein for a discussion of 
the polaron physics in that model) . These energies have been calculated with the pseudo- 
site DMRG method using the corresponding 32 eigenstates as targets. Unfortunately, 
this approach is limited to small number R of targets (of the order of a few tens), which 
is not sufficient for calculating a complete excitation spectrum. 

4T.1. Dynamical correlation functions. The (zero-temperature) linear response of a 
quantum system to a time-dependent perturbation is often given by dynamical correlation 
functions (with h — 1) 

(32) GAiLo + ^v)^--{i'o\A \ ] . -A\^Po), 

TT ho + Lu + irj — H 

where H is the time-independent Hamiltonian of the system, Eq and \ipo) ^-re its ground- 
state energy and wavcfunction, A is the quantum operator corresponding to the physical 
quantity which is analyzed, and is the Hermitian conjugate of A. A small real number 
77 is used in the calculation to shift the poles of the correlation function into the complex 
plane. As a first example, the real part ai{uj > 0) of the optical conductivity is related 
to the imaginary part of the dynamical current-current correlation function, which corre- 
sponds to the operator A = itJ2jai'^]+icr'^j<:r^'^]cr'^j+ia) the above equation for the ID 
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Fig. 10. - Dispersion of the lowest 32 eigenenergies in the one-dimensional Holstein model with 
one electron on a 32-site lattice for ujq = t and g ~ 1 calculated with DMRG (circles). The 
diamonds show a hypothetical polaron dispersion e(fc) = i7p — 2t*[cos(fc) — l], where i?p — —2A71t 
and the effective hopping t* = 0.753f have been fitted to the DMRG dispersion around k = 0. 
The squares show the free-electron dispersion e{k) = — 2tcos(fc). 



Hubbard and Holstein-Hubbard models. As another example, on a one-dimensional lat- 
tice the angle-resolved photocmission spectrum is related to the spectral function A(k, uj) 
given by the imaginary part of with the operator A = J2j e^'^-'cja- 

In general, we are interested in the imaginary part of the correlation function 

(33) Ia{uj + ifi) - Im Ga{oj + *?/) = -{MA^ —— ^— ^— ^A|^o) 

TT [Eo + cj - Hy + if 

in the rj —t Q limit, //i(w) = lim.,,^o Ia{^ + i'n) ^ 0- It should be noted that the spectrum 
Ia{i^ + ivi) for any finite > is equal to the convolution of the spectral function /a(w) 
with a Lorentzian distribution of width rj 

(34) Ia{u + i77) = / dLo'lA{oj')--, L , > 0. 



Let n = 0, 1, 2, . . . be the complete set of eigenstates of H with eigenenergies En 
{\n = 0) corresponds to the ground state \ipo))- The spectral function (|33|) can be written 
in the so-called Lehmann spectral representation 



(35) lA{u; + tv)^-Y.\{n\Am 



[En -Eo- Lof + r?2 



where En — Eq is the excitation energy and |(n|A|0)p the spectral weight of the ra- 
th excited state. Obviously, only states with a finite spectral weight contribute to the 
dynamical correlation function (|32|l and play a role in the dynamics of the physical 
quantity corresponding to the operator A. In the one-dimensional half-filled Hubbard 
model the Hilbert space dimension increases exponentially with the number of sites N 
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but the number of eigenstates with non-zero matrix elements (n|A|0) increases only as a 
power-law of N for the optical conductivity or the one-electron density of states. 

4'1.2. Symmetries . The simplest method for calculating specific excited states uses 
symmetries of the system. If symmetry operators arc well defined in every subsystem, 
DMRG calculations can be carried out to obtain the R lowest eigenstates in a specific 
symmetry subspace. It is also possible to target simultaneously the lowest eigenstates 
in two different subspaces and to calculate matrix elements (n|A|0) between them. This 
allows one to reconstruct the dynamical correlation function at low energy using eq. . 

There are several approaches for using symmetries in DMRG calculations. From a 
computational point of view, the best approach is an explicit implementation of the 
corresponding conserved quantum numbers in the program This is easily done for so- 
called additive quantum numbers such as the particle number or the projection of the 
total spin onto an axis. For instance, if the ^-projection of the total spin is conserved (i.e., 
the spin operator Sz commutes with the Hamilton operator H), one can calculate the 
lowest eigenstates for various quantum numbers Sz to investigate spin excitations |40| . 
As another example, if we study a system with iVc electrons, we can compute the ground- 
state energy Eo{N') for different number A^' of electrons around iVo and thus obtain the 
(charge) gap Egi = Eo{Nc + 1) -f Eo{Nc — 1) — 2Eo{Nc) in the spectrum of free electronic 
charge excitations (see the separate paper by Fehske and Jeckelmann for an application 
of this approach). The extension of DMRG to non-abelian symmetries such as the SU{2) 
spin symmetry of the Hubbard model is presented in ref . |41| . 

A second approach for using symmetries is the construction of projection matrices 
onto invariant subspaces of the Hamiltonian. They can be used to project the matrix 
representation of the superblock Hamiltonian |42) or the initial wave function of the iter- 
ative algorithm used to diagonalizc the superblock Hamiltonian |43| (so that it converges 
to eigenstates in the chosen subspace). This approach has successfully been used to study 
optical excitations in one-dimensional models of conjugated polymers |43[ I44| . 

A third approach consists in adding an interaction term to the system Hamiltonian H 
to shift the eigenstates with the chosen symmetry to lower energies. For instance, if the 
total spin commutes with the Hamilton operator H, one applies the DMRG method 
to the Hamiltonian H' = H + AS^ with A > to obtain the lowest singlet eigenstates 
without interference from the > eigenstates |45| . 

Using symmetries is the most efficient and accurate approach for calculating specific 
low-lying excited states with DMRG. However, its application is obviously restricted 
to those problems which have relevant symmetries and it provides only the lowest R 
eigenstates for given symmetries, where R is at most a few tens for realistic applica- 
tions. Thus this approach cannot describe high-energy excitations nor any complex or 
continuous dynamical spectrum. 

4'1.3. Lanczos-DMRG. The Lanczos-DMRG approach was introduced by Karen 
Hallberg in 1995 as a method for studying dynamical properties of lattice quantum 
many-body systems |46| . It combines DMRG with the continuous fraction expansion 
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technique, also called Lanczos algorithm (see subsubsect. I2"4. l]! . to compute the dynam- 
ical correlation function H32(l . Firstly, the Lanczos algorithm is used to calculate the 
complete dynamical spectrum of a superblock Hamiltonian. Secondly, some Lanczos 
vectors are used as DMRG target in an attempt at constructing an effective Hamilto- 
nian which describes excited states contributing to the correlation function accurately. 
Theoretically, one can systematically improve the accuracy using an increasingly large 
enough number L of Lanczos vectors as targets. Unfortunately, the method becomes nu- 
merically instable for large L as soon as the DMRG truncation error is finite. Therefore, 
Lanczos-DMRG is not an exact numerical method except for a few special cases. The 
cause of the numerical instability is the tendency of the Lanczos algorithm to blow up 
the DMRG truncation errors. In practice, only the first few Lanczos vectors are included 
as target. The accuracy of this type of calculation is unknown and can be very poor. 

From a physical point of view, the failure of the Lanczos-DMRG approach for complex 
dynamical spectra can be understood. With this method one attempts to construct a 
single effective representation of the Hamiltonian H which describes the relevant excited 
states for all excitation energies. This contradicts the essence of a RG calculation, which 
is the construction of an effective representation of a system at a specific energy scale by 
integrating out the other energy scales. 

Nevertheless, Lanczos-DMRG is a relatively simple and quick method for calculating 
dynamical properties within a DMRG approach and it has already been used in several 
works It gives accurate results for systems slightly larger than those which can be 
investigated with exact diagonalization techniques. It is also reliable for larger systems 
with simple discrete spectra made of a few peaks. Moreover, Lanczos-DMRG gives 
accurate results for the first few moments of a spectral function and thus provides us 
with a simple independent check of the spectral functions calculated with other methods. 

In the context of EP systems the Lanczos algorithm has been successfully combined 
with the optimal phonon basis DMRG method (see subsubsect . 13 '4 . 2|l . The optical con- 
ductivity, single-electron spectral functions and electron-pair spectral functions have been 
calculated for the ID Holstein-Hubbard model at various band fillings |^07j. Results 
for these dynamical quantities agree qualitatively with results obtained by conventional 
exact diagonalizations using powerful parallel computers fsect. 

4T.4. Correction vector DMRG. Using correction vectors to calculate dynamical 
correlation functions with DMRG was first proposed by Ramasesha et al. ^48;. The 
correction vector associated with Ga{^ + iv) is defined by 



where |^) = ^I'/'o) is the first Lanczos vector. If the correction vector is known, the 
dynamical correlation function can be calculated directly 



(36) 



|'0A(w + i?7)) 



1 



Eq + uj + irj — H 



(37) 



Ga{oj + if]) 



-(A|V'a(w + *??))■ 
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To calculate a correction vector, one first solves an inhomogeneous linear equation 

(38) [iEo+u;-Hf+rj^]\^)^-rU), 

which always has a unique solution \ip) = \Ya{uj + irf)) for ?7 7^ 0. The correction vector 
is then given by \iPa{'^ + = + ^v)) + + with 

(39) \Xa{uj + iv)) = ^ ' ^° ' \Ya{lu + iJ])). 

V 

One should note that the states \XA{uJ+ir])) and \YA{LO + iT])) are complex if the state \A) 
is not real, but they always determine the real part and imaginary part of the dynamical 
correlation function Ga(w + irj), respectively. The approach can be extended to higher- 
order dynamic response functions such as third-order optical polarizabilities |48j and to 
derivatives of dynamical correlation functions j39) . 

The distinct characteristic of a correction vector approach is that a specific quantum 
state is constructed to compute the dynamical correlation function (|32|) at each 
frequency uj. To obtain a complete dynamical spectrum, the procedure has to be repeated 
for many different frequencies. Therefore, this approach is generally less efficient than 
the iterative methods presented in sect. l2'4l in the context of exact diagonalizations. For 
DMRG calculations, however, this is a highly favorable characteristic. The dynamical 
correlation function can be determined for each frequency to separately using effective 
representations of the system Hamiltonian H and operator A which have to describe a 
single energy scale accurately. 

Kiihner and White 0^1 have devised a correction vector DMRG method which uses 
this characteristic to perform accurate calculations of spectral functions for all frequen- 
cies in large lattice quantum many-body systems. In their method, two correction vectors 
with close frequencies uji and lu2 and finite broadening rj ^ L02 — loi > are included as 
target. The spectrum is then calculated in the frequency interval loi < lu < lu2 using 
eq. (I37|l or the continuous fraction expansion. The calculation is repeated for several 
(overlapping) intervals to determine the spectral function over a large frequency range. 
This procedure makes the accurate computation of complex or continuous spectra pos- 
sible. Nevertheless, there have been relatively few applications of the correction vector 
DMRG method ^Hl because it requires substantial computational resources and is diffi- 
cult to use efficiently. 

4'2. Dynamical DMRG method. - The capability of the correction vector DMRG 
method to calculate continuous spectra shows that using specific target states for each 
frequency is the right approach. Nevertheless, it is highly desirable to simplify this 
approach and to improve its performance. A closer analysis shows that the complete 
problem of calculating dynamical properties can be formulated as a minimization prob- 
lem. This leads to the definition of a more efficient and simpler method, the dynamical 
DMRG (DDMRG) 
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DDMRG enables accurate calculations of dynamical properties for all frequencies in 
large systems with up to a few hundred particles using a workstation. Combined with 
a proper finitc-sizc-scaling analysis (subsect. I4'3|l it also enables the investigation of 
spectral functions in the thermodynamic limit. Therefore, DDMRG provides a powerful 
new approach for investigating the dynamical properties of quantum many-body systems. 

4'2.1. Variational principle. In the correction vector DMRG method the most time- 
consuming task is the calculation of correction vectors in the superblock from cq. (|38|l . 
A well-established approach for solving an inhomogencous linear equation H38|l is to 
formulate it as a minimization problem. Consider the equation Mx = a, where M is 
a positive definite symmetric matrix with a non-degenerate lowest eigenvalue, a is a 
known vector, and x is the unknown vector to be calculated. One can define the function 
M^(x) = X • A/x — X • a — a • X, which has a non-degenerate minimum for the vector 
Xmin = M~^a which is solution of the inhomogencous linear equation. (Kiihncr and 
White |49j used a conjugate gradient method to solve this minimization problem.) 

Generalizing this idea one can formulate a variational principle for dynamical corre- 
lation functions. One considers the functional 

(40) WA,niL^,4') = i^'HEo +uj-Hf+ Tf\i,) + r,{A\i,) + 77(^|A). 

For any rj ^ and a fixed frequency u> this functional has a well-defined and non- 
degenerate minimum for the quantum state which is solution of eq. H38|l . i.e. |'0min) ~ 
\YA{to + irj)) . It is easy to show that the value of the minimum is related to the imaginary 
part of the dynamical correlation function 

(41) WA,T,iui,1pmin) = -TTTjlAiuJ + ilfj . 

Therefore, the calculation of spectral functions can be formulated as a minimization 
problem. To determine lA{'^+i'n) s-t any frequency uj and for any 77 > 0, one minimizes the 
corresponding functional WA,ri{'^,''P)- Once this minimization has been carried out, the 
real part of the correlation function G A{'-^+i'n) can be calculated using eqs. (|37|l and (|39|l . 
This is the variational principle for dynamical correlation functions. It is clear that if we 
can calculate \ YA{^ + iii)) exactly, this variational formulation is completely equivalent to 
a correction vector approach. However, if we can only calculate an approximate solution 
with an error of the order e <C 1, = \Ya{ijJ -\- irf)) + t\4i) with {4>\4>) = 1, the variational 
formulation is more accurate. In the correction vector method the error in the spectrum 
Ia{^ + i'rf) calculated with eq. (|37() is also of the order of e. In the variational approach 
it is easy to show that the error in the value of the minimum WA.ri{^i^mXn)i and thus 
in Ia{'^ + *??), is of the order of e^. With both methods the error in the real part of 
Ga(w -I- irj) is of the order of e. 

The DMRG procedure used to minimize the energy functional E{il]) (see sect.lsj can 
also be used to minimize the functional MOi,r)(^,'0) and thus to calculate the dynamical 
correlation function Ga{^ + iif)- This approach is called the dynamical DMRG method. 
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In principle, it is equivalent to the correction vector DMRG method of Klihner and 
White because the same target states are used to build the DMRG basis in both 
methods. In practice, however, DDMRG has the significant advantage over the correction 
vector DMRG that errors in Iji{u! + iri) are significantly smaller (of the order of instead 
of e as explained above) as soon as DMRG truncation errors are no longer negligible. 

4'2.2. DDMRG algorithm. The minimization of the functional Wa,7j{oj,iP) is easily 
integrated into the usual DMRG algorithm. At every step of a DMRG sweep through the 
system lattice, a superblock representing the system is built and the following calculations 
are performed in the the superblock subspace: 

1. The energy functional Elip) is minimized using a standard iterative algorithm for 
the eigenvalue problem. This yields the ground state vector lipo) a-nd its energy Eq 
in the superblock subspace. 

2. The state |^) is calculated. 

3. The functional WA.riii-^, V') is minimized using an iterative minimization algorithm. 
This gives the first part of the correction vector |yA('j-' + and the imaginary 
part I A [i^ + i??) of the dynamical correlation function through eq. (|41|l . 

4. The second part + irj)) of the correction vector is calculated using eq. (|39|l . 

5. The real part of the dynamical correlation function can be calculated from eq. (|37|1 . 

6. The four states |V'o), IYaI^^ + and + irj)) are included as target in 
the density matrix renormalization to build a new superblock at the next step. 

The robust finite-system DMRG algorithm must be used to perform several sweeps 
through a lattice of fixed size. Sweeps are repeated until the procedure has converged to 
the minimum of both functionals E{ip) and WA,r)(w, V')- 

To obtain the dynamical correlation function GAi^^ + irj) over a range of frequencies, 
one has to repeat this calculation for several frequencies us. If the DDMRG calculations 
are performed independently, the computational effort is roughly proportional to the 
number of frequencies. It is also possible to carry out a DDMRG calculation for several 
frequencies simultaneously, including several states |Xyi(a; + irj)) and |Ya(i^ + ii])) with 
different frequencies oj as target. As calculations for different frequencies are essentially 
independent, it would be easy and very efficient to use parallel computers. 

Because of the variational principle one naively expects that the DDMRG results for 
I A (oj + irj) must converge monotonically from below to the exact result as the number 
m of density matrix eigenstates is increased. In practice, the convergence is less regu- 
lar because of two approximations made to calculate the functional WA^rii^yi^)- First, 
the ground-state energy Eq and the state |^) used in the definition (|40|l of WA,rj{(^,4') 
are not known exactly but calculated with DMRG. Second, one calculates an effective 
representation of H only and assumes that {H'^)eS « (-ffeff)^ to compute WA,rj{t^,i^) in 
the superblock subspace. These approximations can cause a violation of the variational 
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Fig. 11. - Optical conductivity of the one-dimensional half-filled Hubbard model for U — 3t. 
Left panel; DDMRG result calculated on a 128-site lattice using a broadening r/ = O.lt (dashed 
line) and field-theoretical prediction for the infinite chain I38| broadened with a Lorentzian of 
width rj = O.lt (solid line). Right panel: the field-theoretical prediction without broadening 
(solid line) and the DDMRG result after deconvolution (circles). 



bound Wa.tiI'-^ , ip) ^ —i^V^Ai^ + In practice, for a sufficiently large number m of 
density matrix eigenstates kept, the absolute errors in /^(w + irf) decrease systemati- 
cally with increasing m. Errors becomes negligible if enough states are kept to make 
the discarded weight vanish. It is also possible to estimate the accuracy of a DDMRG 
calculation from the results obtained for different values of m. Therefore, DDMRG is an 
exact numerical method as defined in the introdiiction (sect. [T|l . 

The accuracy of the DDMRG approach for ID correlated electron systems has been 
demonstrated by numerous comparisons with exact analytical results |38l 1391 1501 1511 
1521 153L [S^ . As an example, fig. 1111 shows the optical conductivity of the ID half- filled 
Hubbard model for U = 2>t. The DDMRG result (calculated using a broadening 77 = O.lt 
on a 128-site lattice) agrees perfectly with the field-theoretical prediction (also broadened 
with a Lorentzian distribution of width 77) |38l I50| . 

4'3. Spectrum in the thermodynamic limit. - DDMRG allows us to calculate spectral 
functions of a large but finite system with a broadening given by the parameter 77 > 0. To 
determine the properties of a dynamical spectrum in the thermodynamic limit, one has 
to analyze the scaling of the corresponding spectra lN,rii^) as a function of the system 
size N for vanishing broadening rj 

(42) I{lu) = lim lim In,,Aui). 

Computing both limits in this equation from numerical results for /^r j,(ti)) requires a lot of 
accurate data for different values of 77 and N and can be the source of large extrapolation 
errors. A much better approach is to use a broadening r]{N) > which decreases with 
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increasing N and vanishes in the thermodynamic Umit |39| . The dynamical spectrum is 
then given by = Hmjv— >oo lN,r]{N)i^)- From the existence of both Umits in eq. (|42|l it 
can be demonstrated that there exists a minimal broadening %(-/V) > 0, which converges 
to zero for N oo, such that this procedure is exact for all functions ijlN) with 77(A^) > 
770 (TV) and limAr_^oo vi^) = 0- 

The function rio{N) depends naturally on the specific problem studied and can also 
vary for each frequency lu considered. For one-dimensional correlated electron systems 
such as the Hubbard model one finds empirically that a sufficient condition is 

(43) viN) = ^, 

where the constant c is comparable to the effective width of the dynamical spectrum I (to), 
which is finite in such lattice models. This condition has a very simple physical interpre- 
tation. The spectral function I^^^ii^) represents the dynamical response of the system 
over a time period ~ f/77 after one has started to apply an external force. Typically, 
in a lattice model the spectral width is proportional to the velocity of the excitations 
involved in the system response. Thus the condition (|43|l means that excitations are too 
slow to travel the full length ~ of the system in the time interval ^ l/ij and do not 
"sense" that the system is finite. 

An additional benefit of a broadening satisfying the condition H43|) is that the finite- 
system spectrum Ijsi^riiy) becomes indistinguishable from the infinite-system spectrum 
with the same broadening 77 for relatively small N . Therefore, if one knows a spectral 
function I{lo) for an infinite system, its convolution with a Lorentzian of width 7] can 
be compared directly with the numerical results for the finite-system spectrum I^^rii^)- 
This is the approach used in fig. ^] (left panel) to compare DDMRG results for a finite 
lattice with the field-theoretical prediction for an infinite chain. 

Finally, an approximation for an infinite-system (continuous) spectrum can be ob- 
tained by solving the convolution equation (|34|l numerically for /^(w') using the nu- 
merical DDMRG data for a finite system on the left-hand side of this equation |54| . 
Performing such a deconvolution is a ill-conditioned inverse problem, which can only be 
solved approximately using some assumptions on the spectrum properties like its smooth- 
ness. Therefore, the accuracy of deconvolved DDMRG spectra is unknown. In practice, 
however, one obtains often accurate results as shown in fig. 1111 fright panel), where a 
deconvolved DDMRG spectrum is compared to an exact field- theoretical result. 

In summary, the dynamical spectrum of an infinite system can be determined accu- 
rately and efficiently from numerical (DDMRG) data for finite-system spectra using a 
finite-size scaling analysis with a size-dependent broadening ii{N). 

4'4. Application to electron-phonon problems. - The DDMRG algorithm described in 
the previous section can be applied to EP systems such as the Holstein-Hubbard model 
without modification It can naturally be combined with the special DMRG tech- 
niques for systems with bosonic degrees of freedom which are described in subsect. l3"4l As 
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Fig. 12. - Spectral lunctions A{k,uj) (for electron removal, ui < 0, and electron injection, ui > 0) 
of the spinless Holstein model at half filling on a 8-site lattice with periodic boundary conditions. 
The system is in the CDW/Peierls insulating phase (u = O.lt and g — 4). The rapidly oscillating 
thin lines are the KPM results while the smooth thick line are the DDMRG results with the 
pseudo-site method. Note that only jfc| < n/2 is shown because A{k ±n,ij) = A{k, —uj). 



for ground-state simulations, these special techniques substantially reduce the computa- 
tional cost for calculating dynamical properties of EP systems. Nevertheless, applications 
of DDMRG to the Holstein-Hubbard model are significantly more costly than those for 
comparable purely electronic systems such as the Hubbard model. 

As an illustration, we compare DDMRG and ED-KPM (subsubsect. I2"4.2|) results 
for the spectral functions of an eight-site spinless Holstein model in fig. ^| There is an 
excellent overall agreement between both methods. The observable differences are due to 
the larger broadening 77 = O.lt used in the DDMRG simulation. This hides the sharpest 
details of the spectrum like the oscillations which are visible in the KPM spectrum. 

A significant difference between both methods is the computational cost. These 
DDMRG calculations took only 150 CPU hours on an Opteron 244 processor (1.8 GHz) 
and required less than 300 MBytes of memory. It is thus possible to simulate significantly 
larger EP systems than this eight-site lattice with DDMRG while this system size is the 
largest one which can be simulated with exact diagonalization techniques. 

Finally, we note that the broadening 77 = O.lt used for the DDMRG results in fig. 1121 
is one order of magnitude smaller than the value used for a purely electronic systems 
of comparable size (see fig. Ill|l . It seems that the scaling (|43|l is not applicable to EP 
systems. Analytical and numerical investigations will be necessary to determine the 
function rio{N) for EP systems before one can investigate their dynamical properties in 
the thermodynamic limit using the techniques described in subsect. 14 '31 
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5. Conclusion 

Exact diagonalization and density matrix renormalization group techniques are power- 
ful and versatile exact numerical approaches for investigating the properties of electron- 
phonon lattice models for strongly correlated (low-dimensional) materials. Thanks to 
recent developments we can now calculate dynamical quantities which are directly re- 
lated to experimental techniques used in solid-state spectroscopy and often determine 
the properties in the thermodynamic limit using a finite-size scaling analysis. 

* * * 

We would like to thank B. Bauml, H. Benthien, F. Efiler, F. Gebhard, G. Hager, 
S. Nishimoto, G. Wellein, and A. Weisse for valuable discussions. 
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